Björn Reineking, 2024-06-20
Model the relative density of animals (also called range distribution or utilisation distribution) as a function of environmental predictors.
We will use the buffalo data set. Install animove metapackage https://github.com/AniMoveCourse/animove_R_package install.packages(“remotes”) remotes::install_github(“AniMoveCourse/animove_R_package”)
This code needs ctmm version 1.1.1 (or newer) remotes::install_github(“ctmm-initiative/ctmm”) Loading packages
library(animove)
library(ctmm)
library(sf)
library(mvtnorm)
library(terra)
data(buffalo_env)
# Convert from raster to terra object. Terra is the recommended package.
buffalo_env <- rast(buffalo_env)
plot(buffalo_env)
data(buffalo)
projection(buffalo) <- crs(buffalo_env) # project buffalo tracks to projection of raster layers
plot(buffalo_env[[1]])
buffalo_mv <- move2::mt_as_move2(buffalo) # Convert to move2 object for plotting
plot(buffalo_mv, add = TRUE)
## Warning in plot.sf(buffalo_mv, add = TRUE): ignoring all but the first
## attribute
To speed up analyses, we will only work with one individual, Cilla
cilla_mv <- subset(buffalo_mv, track == "Cilla")
plot(buffalo_env[[1]], ext = extent(cilla_mv) * 2)
plot(st_geometry(cilla_mv), add = TRUE, type = "p",
col = "red", pch = 19, cex = 0.5)
Select one individual, Cilla
cilla <- buffalo$Cilla
Fit ctmm model We start with a guesstimate of parameter values; we demand a isotropic ctmm model, because this is a requirement for rsf.fit
cilla_guess <- ctmm.guess(cilla, CTMM = ctmm(isotropic = TRUE),
interactive = FALSE)
# Fit ctmm model and perform model selection, using all but one core
if (file.exists("cilla_fit_rsf.rds")) {
cilla_fit <- readRDS("cilla_fit_rsf.rds")
} else {
cilla_fit <- ctmm.select(cilla, cilla_guess, cores = -1)
saveRDS(cilla_fit, "cilla_fit_rsf.rds")
}
Create named list of rasters rsf.fit can deal with raster layers that differ in resolution or even projection
buffalo_env_list <- list("elev" = raster::raster(buffalo_env[["elev"]]),
"slope" = raster::raster(buffalo_env[["slope"]]),
"var_NDVI" = raster::raster(buffalo_env[["var_NDVI"]]))
Definition of reference grid the akde predictions should get aligned to
reference_grid <- crop(buffalo_env_list[["elev"]], extent(cilla_mv) * 2)
Fit akde
cilla_akde <- akde(cilla, cilla_fit, grid = reference_grid)
## Default grid size of 3 minutes chosen for bandwidth(...,fast=TRUE).
plot(cilla_akde)
# Create a raster representation of the akde, for the reference grid
# PMF is "probability mass function". The values sum up to 1 for the whole grid,
# and represent for a given pixel the probability of finding the animal in that pixel
plot(ctmm::raster(cilla_akde, DF = "PMF"))
The integrator = “Riemann” option is much faster
cilla_rsf_riemann <- rsf.fit(cilla, cilla_akde,
R = buffalo_env_list, integrator = "Riemann")
## Warning in .local(x, ...): This function is only useful for Raster* objects
## with a longitude/latitude coordinates
## Maximizing likelihood @ n=164976
## Calculating Hessian
# Estimates and confidence intervals for all model parameters
# - Selection: The three environmental layers
# - Movement: home range area,time scales of position and velocity autocorrelation;
# derived quantities speed and diffusion
summary(cilla_rsf_riemann)
## $name
## [1] "OUF"
##
## $DOF
## mean area diffusion speed
## 8.60272 17.18711 74.92798 74.34974
##
## $CI
## low est high
## var_NDVI (1/var_NDVI) -186.50622449 -73.08023745 40.34574959
## slope (1/slope) -0.37987332 -0.03789964 0.30407403
## elev (1/elev) -0.04621166 -0.01264287 0.02092592
## area (square kilometers) 232.61241367 397.97433177 607.01231642
## τ[position] (days) 4.47294577 7.59098047 12.88255826
## τ[velocity] (minutes) 40.67461444 43.17064097 45.81983793
## speed (kilometers/day) 12.10304972 13.65429591 15.20333307
## diffusion (square kilometers/day) 4.28602881 5.44972358 6.75106371
A suitability map - with 95% confidence interval
suitability_riemann <- ctmm::suitability(cilla_rsf_riemann, R = buffalo_env_list,
grid = reference_grid)
terra::plot(suitability_riemann)
Range distribution (includes the ranging behaviour)
agde_cilla <- agde(CTMM = cilla_rsf_riemann, R = buffalo_env_list,
grid = reference_grid)
plot(agde_cilla)
# Plot raster of range distribution
agde_raster <- ctmm::raster(agde_cilla, DF = "PMF")
plot(agde_raster)
Selection-informed akde
akde_rsf <- akde(cilla, CTMM = cilla_rsf_riemann, R = buffalo_env_list,
grid = reference_grid)
## Default grid size of 3 minutes chosen for bandwidth(...,fast=TRUE).
plot(ctmm::raster(akde_rsf, DF = "PMF"))
# Plot all three probability mass functions with a common color scale
pmf <- raster::stack(ctmm::raster(cilla_akde, DF = "PMF"),
ctmm::raster(agde_cilla, DF = "PMF"),
ctmm::raster(akde_rsf, DF = "PMF"))
names(pmf) <- c("akde", "rsf.fit", "akde_rsf.fit")
plot(pmf, zlim=c(0, max(raster::getValues(pmf)))) # use zlim to force same color scale
Functions to generate quadrature points and predict with the model
rsf_points <- function(x, UD, R = NULL, n = 1e5, k = 1e6, type = "Riemann",
rmax = 6*sqrt(UD@CTMM$sigma[1,1]),
interpolation = FALSE) {
# Samples background points from a 2D normal distribution fitted to the relocation data, and extracts environmental
# information from a raster object "R"
# x: telemetry object
# UD: UD object
# R: terra object
# n: number of background points to sample
# k: weight of presence points
# rmax: maximum distance for Riemann-type integration
# interpolation: do interpolation when sampling the grid
# When type 0´= "MonteCarlo", importance sampling is done
stopifnot(UD@CTMM$isotropic)
stopifnot(type %in% c("Riemann", "MonteCarlo"))
if (type == "Riemann") {
quadrature_pts <- values(R) #raster::getValues(R)
xy <- as.data.frame(xyFromCell(R, 1:ncell(R)))
xy <- sf::st_as_sf(xy, coords = c("x", "y"), crs = crs(R))
xy <- st_transform(xy, UD@info$projection)
xy <- st_coordinates(xy)
r <- sqrt(((xy[,1] - UD@CTMM$mu[1]))^2 + ((xy[,2] - UD@CTMM$mu[2]))^2)
bg <- data.frame(case_ = 0,
x_ = xy[r<rmax,1], y_ = xy[r<rmax,2], w_ = prod(res(R)), k_ = k)
bg <- cbind(bg, quadrature_pts[r<rmax,])
bg <- sf::st_as_sf(bg, coords = c("x_", "y_"), crs = UD@info$projection)
xx <- data.frame(case_ = 1, x_ = x$x, y_ = x$y,
w_ = 1/k * UD$weights * mean(UD$DOF.area), k_ = k)
xx <- sf::st_as_sf(xx, coords = c("x_", "y_"), crs = UD@info$projection)
xx[names(R)] <- as.data.frame(extract(R, st_transform(xx, crs(R)),
ID = FALSE,
method = ifelse(interpolation, "bilinear", "simple")))
xx <- rbind(bg, xx)
} else {
quadrature_pts <- MASS::mvrnorm(n, mu = UD@CTMM$mu, Sigma = UD@CTMM$sigma)
xx <- data.frame(case_ = 0, x_ = quadrature_pts[, 1], y_ = quadrature_pts[, 2], w_ = UD@CTMM$sigma[1,1]/n, k_ = k)
xx <- rbind(xx, data.frame(case_ = 1, x_ = x$x, y_ = x$y,
w_ = 1/k * UD$weights * mean(UD$DOF.area), k_ = k
))
xx <- sf::st_as_sf(xx, coords = c("x_", "y_"), crs = UD@info$projection)
xx[names(R)] <- as.data.frame(extract(R, st_transform(xx, crs(R)),
ID = FALSE,
method = ifelse(interpolation, "bilinear", "simple")))
}
xy <- st_coordinates(xx)
colnames(xy) <- c("x_", "y_")
sd <- sqrt(UD@CTMM$sigma[1,1])
xy[,1] <- (xy[,1] - UD@CTMM$mu[1])/sd
xy[,2] <- (xy[,2] - UD@CTMM$mu[2])/sd
xx <- cbind(xy, xx)
xx
}
predict_rsf <- function(model, UD, object, include_avail = TRUE, data_crs = UD@info$projection) {
if(include_avail) {
xy <- as.data.frame(xyFromCell(object, 1:ncell(object)))
xy <- sf::st_as_sf(xy, coords = c("x", "y"), crs = crs(object))
xy <- st_transform(xy, data_crs)
xy <- st_coordinates(xy)
colnames(xy) <- c("x_", "y_")
sd <- sqrt(UD@CTMM$sigma[1,1])
xy[,1] <- (xy[,1] - UD@CTMM$mu[1])/sd
xy[,2] <- (xy[,2] - UD@CTMM$mu[2])/sd
} else {
xy <- cbind("x_" = rep(0, ncell(object)), "y_" = rep(0, ncell(object)))
}
newdata <- as.data.frame(cbind(xy, values(object)))
lambda <- as.numeric(predict(model, newdata, type = "link"))
r <- rast(object[[1]])
r[] <- exp(lambda - max(lambda, na.rm = TRUE))
r <- r / sum(values(r), na.rm = TRUE)
r
}
Generate quadrature points (“background points”)
set.seed(2)
rsf_cilla_df <- rsf_points(cilla, cilla_akde, buffalo_env, interpolation = TRUE)
rsf_cilla_df <- rsf_cilla_df[!is.na(rsf_cilla_df$slope),] # Remove lines with NA values
Fit a downweighted Poisson regression the homeranging behaviour is represented by x_ + y_ + I(-(x_^2 + y_^2)/2)
m_rsf_cilla <- glm(case_*k_ ~ x_ + y_ + I(-(x_^2 + y_^2)/2) + elev + slope + var_NDVI,
family = poisson(), data= rsf_cilla_df, weights = w_)
## Warning: glm.fit: fitted rates numerically 0 occurred
Summary of model and confidence intervals for parameter estimates
summary(m_rsf_cilla)
##
## Call:
## glm(formula = case_ * k_ ~ x_ + y_ + I(-(x_^2 + y_^2)/2) + elev +
## slope + var_NDVI, family = poisson(), data = rsf_cilla_df,
## weights = w_)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -10.29200 3.50561 -2.936 0.00333 **
## x_ 0.30505 0.35459 0.860 0.38963
## y_ 0.48068 0.44149 1.089 0.27626
## I(-(x_^2 + y_^2)/2) 1.14760 0.27568 4.163 3.14e-05 ***
## elev -0.01262 0.01711 -0.737 0.46089
## slope -0.03798 0.17444 -0.218 0.82766
## var_NDVI -73.10604 57.87565 -1.263 0.20653
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 1164.6 on 290412 degrees of freedom
## Residual deviance: 1093.8 on 290406 degrees of freedom
## AIC: 1107.8
##
## Number of Fisher Scoring iterations: 24
confint(m_rsf_cilla)
## 2.5 % 97.5 %
## (Intercept) -16.51698320 -3.01359591
## x_ -0.37703096 1.03865127
## y_ -0.32870578 1.39531290
## I(-(x_^2 + y_^2)/2) 0.68634033 1.76974590
## elev -0.04798511 0.01659574
## slope -0.48589830 0.19826987
## var_NDVI -184.25888765 43.40197008
Map of suitability, including the home ranging behaviour
suitability_glm <- predict_rsf(m_rsf_cilla, cilla_akde,
crop(buffalo_env, extent(reference_grid)))
plot(suitability_glm)
Map of suitability, without the home ranging behaviour
suitability_no_avail_glm <- predict_rsf(m_rsf_cilla, cilla_akde,
crop(buffalo_env, extent(reference_grid)),
include_avail = FALSE)
plot(suitability_no_avail_glm)
Comparison of estimates from the two approaches (ctmm and “classic” via glm)
vars <- c("elev", "slope", "var_NDVI")
# Estimates are not identical but similar
cilla_rsf_riemann$beta[vars]
## elev slope var_NDVI
## -0.01264287 -0.03789964 -73.08023745
coef(m_rsf_cilla)[vars]
## elev slope var_NDVI
## -0.01261716 -0.03797706 -73.10604323
# Confidence intervals are of similar width, but location is different
ci_glm <- confint(m_rsf_cilla)
ci_glm[vars, ]
## 2.5 % 97.5 %
## elev -0.04798511 0.01659574
## slope -0.48589830 0.19826987
## var_NDVI -184.25888765 43.40197008
summary(cilla_rsf_riemann)$CI[sapply(vars, function(x) grep(x, rownames(summary(cilla_rsf_riemann)$CI))), ]
## low est high
## elev (1/elev) -0.04621166 -0.01264287 0.02092592
## slope (1/slope) -0.37987332 -0.03789964 0.30407403
## var_NDVI (1/var_NDVI) -186.50622449 -73.08023745 40.34574959
Create maps of suitability based on parameter estimates
M <- terra::values(buffalo_env)
eta_glm <- M[,vars] %*% coef(m_rsf_cilla)[vars]
eta_rsf <- M[,vars] %*% cilla_rsf_riemann$beta[vars]
suit_r_glm <- buffalo_env[[1]]
suit_r_glm[] <- exp(eta_glm)
suit_r_rsf <- buffalo_env[[1]]
suit_r_rsf[] <- exp(eta_rsf)
suit_raster <- c(suit_r_rsf, suit_r_glm)
names(suit_raster) <- c("ctmm", "classic")
plot(suit_raster)
plot(suit_raster, ext = extent(reference_grid))
h_avail: Histogram of elevation values within the reference grid (here taken as available elevations) h_use: Histogram of elevation values at locations used by Cilla
h_avail <- hist(values(reference_grid), prob = TRUE, col=rgb(0,0,1,1/4), main = "",
xlab = "Elevation (m a.s.l.)", ylim = c(0,0.015))
h_use <- hist(extract(reference_grid, cilla_mv), prob = TRUE,
breaks = h_avail$breaks, add = TRUE, col=rgb(1,0,0,1/4))
legend("topright",
fill = c(rgb(0,0,1,1/4), rgb(1,0,0,1/4)),
legend = c("available","use"),
bty = "n"
)
# A simple rsf model containing only selection for elevation, but with
# a quadratic relationship. No homeranging is included
m_rsf_cilla_elev <- glm(case_*k_ ~ poly(elev, 2), family = poisson(),
data= rsf_cilla_df, weights = w_)
## Warning: glm.fit: fitted rates numerically 0 occurred
# Estimate of density of available elevation.
# We need this to calculate the integration constant to convert the
# habitat selection values (i.e. the prediction of the habitat model)
# to the use/availability ratio.
d_avail <- density(values(reference_grid))
elev_pred <- data.frame("elev" = d_avail$x)
elev_pred$w <- predict(m_rsf_cilla_elev,
newdata = elev_pred,
type = "response")
Calculate the integration constant \[u(x) = \frac{w(x) a(x)}{\int_{E}w(X) a(X)dX} \] Where u(x): use of resource level x ; a(x): availablility of resource level x w(x): selection of resource level x, E: the range of resource levels. K is the integration constant \(\int_{E}w(X) a(X)dX\).
K <- sum(elev_pred$w * d_avail$y * diff(d_avail$x[1:2]))
plot(h_avail$mids, h_use$density/h_avail$density,
xlab = "Elevation (m a.s.l.)",
ylab = "use/availability ratio")
abline(h = 1, lty = 3)
lines(elev_pred$elev, elev_pred$w / K, col = "green", lwd = 2)
legend("topright",
col = c("black", "green"),
lty = c(NA,1), lwd = 2, pch = c(1, NA),
legend = c("use/availability","rsf: poly(elev, 2)/K"),
bty = "n"
)